library(cluster)
library(broom)
library(glmnet)
library(modelr)
library(ggplot2)
library(cowplot)
library(mongolite)
library(ggmap)
library(dplyr)
library(sp)
library(lubridate)
library(tidyr)
library(reshape2)
library(stringr)
preciosModelo <- read.csv(file = '/home/ignacio/datos/facultad/repos/tpEspecializacion/data/preciosModelo.csv')
preciosModelo = select (preciosModelo,-c(X,barrio,sucursal))
# banderaDescripcion + medicion + barrio + banderaDescripcion + pVentasC + mCuadradoC
lm_precio2bandera = lm(formula = precio~banderaDescripcion, data=preciosModelo)
lm_precio2medicion = lm(formula = precio~medicion, data=preciosModelo)
#lm_precio2barrio = lm(formula = precio~barrio, data=precios)
summary(lm_precio2bandera)
##
## Call:
## lm(formula = precio ~ banderaDescripcion, data = preciosModelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82.83 -38.71 -18.17 14.82 617.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 81.1674 0.1042 778.923 < 2e-16
## banderaDescripcionDisco 5.0769 0.1752 28.973 < 2e-16
## banderaDescripcionExpress -5.3563 0.2768 -19.354 < 2e-16
## banderaDescripcionHipermercado Carrefour -3.8789 0.3155 -12.294 < 2e-16
## banderaDescripcionJOSIMAR SUPERMERCADOS -5.3012 0.8146 -6.507 7.65e-11
## banderaDescripcionJumbo 4.7310 0.3405 13.894 < 2e-16
## banderaDescripcionMarket -3.9604 0.1744 -22.714 < 2e-16
## banderaDescripcionMi Changomas -10.5035 1.0131 -10.368 < 2e-16
## banderaDescripcionSupermercados DIA -4.6258 0.2690 -17.198 < 2e-16
## banderaDescripcionVea 1.9413 0.2553 7.605 2.85e-14
##
## (Intercept) ***
## banderaDescripcionDisco ***
## banderaDescripcionExpress ***
## banderaDescripcionHipermercado Carrefour ***
## banderaDescripcionJOSIMAR SUPERMERCADOS ***
## banderaDescripcionJumbo ***
## banderaDescripcionMarket ***
## banderaDescripcionMi Changomas ***
## banderaDescripcionSupermercados DIA ***
## banderaDescripcionVea ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 65.55 on 1146544 degrees of freedom
## Multiple R-squared: 0.003011, Adjusted R-squared: 0.003003
## F-statistic: 384.7 on 9 and 1146544 DF, p-value: < 2.2e-16
#coef(lm_precio2bandera)
summary(lm_precio2medicion)
##
## Call:
## lm(formula = precio ~ medicion, data = preciosModelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -81.69 -38.51 -18.25 15.34 610.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.23229 0.13200 569.95 <2e-16 ***
## medicion 1.01898 0.02121 48.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 65.58 on 1146552 degrees of freedom
## Multiple R-squared: 0.002008, Adjusted R-squared: 0.002007
## F-statistic: 2307 on 1 and 1146552 DF, p-value: < 2.2e-16
glance(lm_precio2bandera)
glance(lm_precio2medicion)
# banderaDescripcion + sucursalTipo + medicion + pVentasC + mCuadradoC
lm_precioMultiple = lm(precio ~ banderaDescripcion + sucursalTipo + medicion + pVentasC + mCuadradoC , data=preciosModelo)
# medicion - barrio - sucursalTipo - banderaDescripcion
summary(lm_precioMultiple)
##
## Call:
## lm(formula = precio ~ banderaDescripcion + sucursalTipo + medicion +
## pVentasC + mCuadradoC, data = preciosModelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -87.35 -38.54 -17.93 15.02 615.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.89405 0.20196 375.787 < 2e-16
## banderaDescripcionDisco 5.38072 0.20149 26.705 < 2e-16
## banderaDescripcionExpress -4.89703 0.29251 -16.741 < 2e-16
## banderaDescripcionHipermercado Carrefour -4.05007 0.34297 -11.809 < 2e-16
## banderaDescripcionJOSIMAR SUPERMERCADOS -4.87521 0.83306 -5.852 4.85e-09
## banderaDescripcionJumbo 4.89757 0.34393 14.240 < 2e-16
## banderaDescripcionMarket -3.77411 0.18581 -20.312 < 2e-16
## banderaDescripcionMi Changomas -9.39337 1.02765 -9.141 < 2e-16
## banderaDescripcionSupermercados DIA -4.28939 0.28543 -15.028 < 2e-16
## banderaDescripcionVea 2.29485 0.26545 8.645 < 2e-16
## sucursalTipoSupermercado -0.50046 0.17477 -2.864 0.004188
## medicion 1.01813 0.02119 48.052 < 2e-16
## pVentasCbajo 0.34679 0.26761 1.296 0.195011
## pVentasCmedio 0.12908 0.17839 0.724 0.469339
## mCuadradoCbajo 0.02740 0.18938 0.145 0.884940
## mCuadradoCmedio -0.52743 0.14472 -3.645 0.000268
##
## (Intercept) ***
## banderaDescripcionDisco ***
## banderaDescripcionExpress ***
## banderaDescripcionHipermercado Carrefour ***
## banderaDescripcionJOSIMAR SUPERMERCADOS ***
## banderaDescripcionJumbo ***
## banderaDescripcionMarket ***
## banderaDescripcionMi Changomas ***
## banderaDescripcionSupermercados DIA ***
## banderaDescripcionVea ***
## sucursalTipoSupermercado **
## medicion ***
## pVentasCbajo
## pVentasCmedio
## mCuadradoCbajo
## mCuadradoCmedio ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 65.48 on 1146538 degrees of freedom
## Multiple R-squared: 0.005041, Adjusted R-squared: 0.005028
## F-statistic: 387.2 on 15 and 1146538 DF, p-value: < 2.2e-16
precioMultiple_resid = augment(lm_precioMultiple)
precioMultiple_resid
#El promedio de los residuos debe ser un numero muy cercano a cero
mean(precioMultiple_resid$.resid)
## [1] 1.106236e-11
Como se puede apreciar el valor obtenido del promedio de todos los residuos, es un numero cercano a cero.
ggplot(precioMultiple_resid, aes(precioMultiple_resid$.resid)) +
geom_freqpoly(binwidth = 1.5)+
labs(fill = "precioMultiple_resid$.resid", title = "Poligono de frecuencia de los residuos", x = "Residuo", y = "count")
## Warning: Use of `precioMultiple_resid$.resid` is discouraged. Use `.resid`
## instead.
ggplot(precioMultiple_resid, aes(sample= .std.resid))+
stat_qq()+
geom_abline(colour = "blue", size=2)+
labs(title = "QQ plot Normal(0,1)",subtitle = "Regresion lineal multiple", x = "Valores teóricos", y = "Residuos estandarizados")+
theme_classic(base_size = 25)
Se quiere validar, si los residuos siguien una distribucion teorica, N(0,1). Como podemos ver el modelo en los extremos tiende a alejarse de la distribucion Normal, por lo que puedo concluir que el modelo no esta bien definido.
ggplot(precioMultiple_resid, aes(.fitted, .resid)) +
geom_point()+
geom_hline(yintercept = 0) +
geom_smooth(se = FALSE)+
labs(title = "Residuos versus el modelo ajustado",subtitle = "Regresion lineal multiple", x = "valores fitted", y = "Residuos")+
theme_classic(base_size = 20)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Otro caso interesando para estudiar, es si los residuos tienen o no una estructura definida. Lo que se obseva es una clara estructura en el medio del grafico, esto esta indicando que una parte sistemÔtica del fenómeno que se esta perdiendo, lo cual indica que el modelo no esta funcionando como se esperaria.
# banderaDescripcion + sucursalTipo + medicion + banderaDescripcion + pVentasC + mCuadradoC
preciosModelo_log = preciosModelo
preciosModelo_log$precio = log(preciosModelo_log$precio)
preciosModelo_log$medicion = log(preciosModelo_log$medicion)
lm_precioMultiple_log = lm(precio ~ banderaDescripcion + sucursalTipo + medicion + pVentasC + mCuadradoC , data=preciosModelo_log)
summary(lm_precioMultiple_log)
##
## Call:
## lm(formula = precio ~ banderaDescripcion + sucursalTipo + medicion +
## pVentasC + mCuadradoC, data = preciosModelo_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.09431 -0.40124 -0.00623 0.41589 2.43126
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 4.0812515 0.0022879 1783.825
## banderaDescripcionDisco 0.0615033 0.0021787 28.230
## banderaDescripcionExpress -0.0516603 0.0031628 -16.334
## banderaDescripcionHipermercado Carrefour -0.0489795 0.0037083 -13.208
## banderaDescripcionJOSIMAR SUPERMERCADOS -0.0282972 0.0090073 -3.142
## banderaDescripcionJumbo 0.0580706 0.0037187 15.616
## banderaDescripcionMarket -0.0463983 0.0020091 -23.095
## banderaDescripcionMi Changomas -0.0691800 0.0111116 -6.226
## banderaDescripcionSupermercados DIA -0.0381283 0.0030862 -12.354
## banderaDescripcionVea 0.0280666 0.0028702 9.778
## sucursalTipoSupermercado -0.0055061 0.0018896 -2.914
## medicion 0.0473387 0.0009461 50.035
## pVentasCbajo 0.0033629 0.0028935 1.162
## pVentasCmedio 0.0005203 0.0019288 0.270
## mCuadradoCbajo -0.0009798 0.0020476 -0.479
## mCuadradoCmedio -0.0064328 0.0015648 -4.111
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## banderaDescripcionDisco < 2e-16 ***
## banderaDescripcionExpress < 2e-16 ***
## banderaDescripcionHipermercado Carrefour < 2e-16 ***
## banderaDescripcionJOSIMAR SUPERMERCADOS 0.00168 **
## banderaDescripcionJumbo < 2e-16 ***
## banderaDescripcionMarket < 2e-16 ***
## banderaDescripcionMi Changomas 4.79e-10 ***
## banderaDescripcionSupermercados DIA < 2e-16 ***
## banderaDescripcionVea < 2e-16 ***
## sucursalTipoSupermercado 0.00357 **
## medicion < 2e-16 ***
## pVentasCbajo 0.24514
## pVentasCmedio 0.78735
## mCuadradoCbajo 0.63229
## mCuadradoCmedio 3.94e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.708 on 1146538 degrees of freedom
## Multiple R-squared: 0.00547, Adjusted R-squared: 0.005457
## F-statistic: 420.4 on 15 and 1146538 DF, p-value: < 2.2e-16
lm_precioMultiple_log_resid = augment(lm_precioMultiple_log)
lm_precioMultiple_log_resid
mean(lm_precioMultiple_log_resid$.resid)
## [1] 5.4611e-12
ggplot(lm_precioMultiple_log_resid, aes(lm_precioMultiple_log_resid$.resid)) +
geom_freqpoly(binwidth = 2.5)+
labs(fill = "propiedades_resid$.resid", title = "Poligono de frecuencia de los residuos", x = "Residuo", y = "count")
## Warning: Use of `lm_precioMultiple_log_resid$.resid` is discouraged. Use
## `.resid` instead.
ggplot(lm_precioMultiple_log_resid, aes(sample= .std.resid))+
stat_qq()+
geom_abline(colour = "blue", size=2)+
labs(title = "QQ plot Normal(0,1)",subtitle = "Regresion lineal multiple logaritmica", x = "Valores teóricos", y = "Residuos estandarizados")+
theme_classic(base_size = 25)
Lo que se obsera en este grafico, es que si bien en los extremos la tendencia es alejarse de la recta, los valores estan mucho mas pegados a ella que en el modelo anterior, lo mismo ocurre con los valores intermedios que estan practicamente sobre la recta. Por lo antes explicado, este modelo esta mejor definido que el anterior.
ggplot(lm_precioMultiple_log_resid, aes(.fitted, .resid)) +
geom_point()+
geom_hline(yintercept = 0) +
geom_smooth(se = FALSE)+
labs(title = "Residuos versus el modelo ajustado",subtitle = "Regresion lineal multiple logaritmica", x = "valores fitted", y = "Residuos")+
theme_classic(base_size = 20)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Si bien en este caso la diferencia no es tan notoria como en el analisis anterior, se puede apreciar que los residuos no estan formando una figura tan concentrada con en el caso no logaritmico, dando una mejora al modelo en este caso. Repasando el articulo sobre la aplicacion de logaritmos para el estudio, este nuevo modelo con logaritmos podria considerarse un hibrido entre un modelo log-nivel para las covariables que no se modificaron y un modelo log-log para aquellas que si lo fueron.
lineal_coef= lm_precioMultiple %>% tidy(conf.int=TRUE)
lineal_coef_log= lm_precioMultiple_log %>% tidy(conf.int=TRUE)
ggplot(lineal_coef, aes(term, estimate))+
geom_point()+
geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
labs(title = "Coeficientes de la regresion lineal", x="", y="Estimacion e Int. Confianza") +
theme_bw() +
theme(axis.text.x = element_text(angle=90))
ggplot(lineal_coef_log, aes(term, estimate))+
geom_point()+
geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
labs(title = "Coeficientes de la regresion lineal log", x="", y="Estimacion e Int. Confianza") +
theme_bw() +
theme(axis.text.x = element_text(angle=90))
ggplot(lineal_coef, aes(reorder(term, -p.value), p.value, fill=p.value))+
geom_bar(stat = 'identity', aes(fill=p.value))+
geom_hline(yintercept = 0.05) +
labs(title = "P-valor de los regresores para multiple", x="", y="P-valor") +
theme_bw() +
theme(axis.text.x = element_text(angle=90)) +
scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )
ggplot(lineal_coef_log, aes(reorder(term, -p.value), p.value, fill=p.value))+
geom_bar(stat = 'identity', aes(fill=p.value))+
geom_hline(yintercept = 0.05) +
labs(title = "P-valor de los regresores para multiple log", x="", y="P-valor") +
theme_bw() +
theme(axis.text.x = element_text(angle=90)) +
scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )
multiple = lm_precioMultiple %>% glance() %>% select(r.squared, adj.r.squared, p.value)
multiple_log = lm_precioMultiple_log %>% glance() %>% select(r.squared, adj.r.squared, p.value)
bind_rows(multiple, multiple_log) %>% mutate(modelo= c('multiple', 'multiple_log'))
train_test <- preciosModelo %>% resample_partition(c(train=0.7,test=0.3))
precios_train <- train_test$train %>% as_tibble()
precios_test <- train_test$test %>% as_tibble()
# Vector con los salarios
prod_precios = preciosModelo$precio
# Matriz con los regresores
prod_mtx = model.matrix(precio~ banderaDescripcion + sucursalTipo + medicion + pVentasC + mCuadradoC, data = preciosModelo)
# Modelo Lasso
lasso.mod=glmnet(x=prod_mtx, # Matriz de regresores
y=prod_precios, #Vector de la variable a predecir
alpha=1, # Indicador del tipo de regularizacion
standardize = F) # Que esta haciendo este parametro?
lasso_coef = lasso.mod %>% tidy()
lasso_coef
plot(lasso.mod, xvar = "lambda", label = TRUE)
title(main = list("Evolución de los coeficientes Lasso",
cex = 1.5, col = "black", font = 10))
#plot(lasso.mod)
Evolución de los coeficientes a medida que se incrementa λ. Los coeficientes se van haciendo mÔs pequeños a medida que se incrementa el valor de λ.
# Graficos para los valores de lambda en ggplot.
g1=lasso_coef %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Lasso con Intercepto", y="Coeficientes")
g2=lasso_coef %>% filter(term!='(Intercept)') %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Lasso sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)
lasso_cv=cv.glmnet(x=prod_mtx,y=prod_precios,alpha=1, standardize = T)
lasso_cv
##
## Call: cv.glmnet(x = prod_mtx, y = prod_precios, alpha = 1, standardize = T)
##
## Measure: Mean-Squared Error
##
## Lambda Measure SE Nonzero
## min 0.0212 4288 13.01 14
## 1se 1.8476 4301 12.82 2
plot(lasso_cv)
title(main = list("Lasso",
cex = 2.5, col = "black", font = 10))
El grƔfico nos muestra la media del MSE con su limite superior e inferior y la cantidad de varaibles que sobreviven para cada valor de lambda.
# Información de CV en dataframe con tidy
lasso_cv %>% tidy()
# Lambda minimo y lambda a 1 desvio estandar
lasso_cv %>% glance()
# Selección lambda óptimo
lasso_lambda_opt = lasso_cv$lambda.min
# Entrenamiento modelo óptimo
lasso_opt = glmnet(x=prod_mtx, # Matriz de regresores
y=prod_precios, #Vector de la variable a predecir
alpha=1, # Indicador del tipo de regularizacion
standardize = TRUE, # Estandarizamos
lambda = lasso_lambda_opt)
# Salida estandar
#lasso_opt
# Tidy
lasso_opt %>% tidy()
Hay quedado 14 variables explicando el 0.5134 % del deviance.
#Modelo ridge
ridge.mod=glmnet(x=prod_mtx, # Matriz de regresores
y=prod_precios, #Vector de la variable a predecir
alpha=0, # Indicador del tipo de regularizacion
standardize = TRUE)
#Coeficientes tidy
ridge_coef= ridge.mod %>% tidy()
ridge_coef
#plot(ridge.mod)
plot(ridge.mod, xvar = "lambda", label = TRUE)
title(main = list("Evolución de los coeficientes Ridge",
cex = 1.5, col = "black", font = 10))
g1=ridge_coef %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Ridge con Intercepto", y="Coeficientes")
g2=ridge_coef %>% filter(term!='(Intercept)') %>%
ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw() + theme(legend.position = 'none') +
labs(title="Ridge sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)
ridge_cv=cv.glmnet(x=prod_mtx,y=prod_precios,alpha=0, standardize = T)
plot(ridge_cv)+
title(main = list("Ridge",
cex = 2.5, col = "black", font = 10))
## integer(0)
# Selección lambda óptimo
lasso_cv %>% glance()
# Selección lambda óptimo
ridge_lambda_opt = ridge_cv$lambda.min
# Entrenamiento modelo óptimo
ridge_opt = glmnet(x=prod_mtx, # Matriz de regresores
y=prod_precios, #Vector de la variable a predecir
alpha=0, # Indicador del tipo de regularizacion
standardize = TRUE, # Estandarizamos
lambda = ridge_lambda_opt)
# Salida estandar
#ridge_opt
ridge_opt %>% tidy()
ridge_dev = ridge_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
ggplot(., aes(log(lambda), dev.ratio)) +
geom_point() +
geom_line() +
geom_vline(xintercept = log(ridge_lambda_opt), color='steelblue', size=1.5) +
labs(title='Ridge: Deviance') +
theme_bw()
lasso_dev = lasso_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
ggplot(., aes(log(lambda), dev.ratio)) +
geom_point() +
geom_line() +
geom_vline(xintercept = log(lasso_lambda_opt), color='firebrick', size=1.5) +
labs(title='Lasso: Deviance') +
theme_bw()
plot_grid(ridge_dev, lasso_dev)
Compracion de la relación entre el porcentaje de deviance explicada y lambda para los tres tipos de modelos que realizamos
ridge_opt
##
## Call: glmnet(x = prod_mtx, y = prod_precios, alpha = 0, lambda = ridge_lambda_opt, standardize = TRUE)
##
## Df %Dev Lambda
## 1 15 0.5 0.4684
eval_results <- function(true, predicted, df) {
SSE <- sum((predicted - true)^2)
SST <- sum((true - mean(true))^2)
R_square <- 1 - SSE / SST
RMSE = sqrt(SSE/nrow(df))
# Model performance metrics
data.frame(
RMSE = RMSE,
Rsquare = R_square
)
}
# Prediccion y evaluacion en train data Lasso
predictions_train <- predict(lasso_opt, s = lasso_lambda_opt, newx = prod_mtx)
eval_results(preciosModelo$precio, predictions_train, preciosModelo)
# Prediccion y evaluacion en test data Lasso
predictions_test <- predict(lasso_opt, s = lasso_lambda_opt, newx = prod_mtx)
eval_results(preciosModelo$precio, predictions_test, preciosModelo)
# Prediction and evaluation on train data Ridge
predictions_train <- predict(ridge_opt, s = ridge_lambda_opt, newx = prod_mtx)
eval_results(preciosModelo$precio, predictions_train, preciosModelo)
# Prediccion y evaluacion en test data Ridge
predictions_test <- predict(ridge_opt, s = ridge_lambda_opt, newx = prod_mtx)
eval_results(preciosModelo$precio, predictions_test, preciosModelo)